GEMC simulation results for custom Helium potential

Eliseo,

Thank you for responding to me question. I need to apologize though, I thought that Ed had removed this question after I emailed on August 22nd saying that I found this energy units issue. Indeed, by converting to amu A^2/ps^2 we get the correct GEMC results. I hope that you did not spend too much time trying to debug our code.

Thanks again

Eliseo and/or Ed,

Now that we have the two-body Helium potential working, our next (and final) task is to implement a three-body Helium potential.

We already have the Fortran 90 code written that takes as input three interatomic distances (rij, rik, and rjk) and outputs the three-body energy. It was readily obvious for the two-body Helium potential that we just needed to call our two-body function in the Compute_AtomPair_Energy subroutine. However, since Cassandra does not have a "Compute_ThreeBody_Energy" subroutine, it is not as clear to me what the best approach is moving forward. Since you are much more familiar with the Cassandra code and have likely given some thought to the implementation of three-body potentials, I was hoping you could provide some insight.

Here is my (naive) idea so far:

1) After calling Compute_MolecularPair_Energy within Compute_Molecule_Nonbond_Inter_Energy (around line 884 in energy_routines), include an additional loop over all the molecules/atoms that are not "is" or "ispecies" (call this "ithird")
2) Call Check_MolecularPair_Cutoff for the distances between ithird and is/ispecies.
3) If these additional two distances are also within the cut-off, call Compute_ThreeBody_Energy and add this to the Eij_inter_vdw term (we could try separating this into its own variable if we want to)

Note that I'm not sure right now how we would avoid double counting the interactions with this approach.

One benefit of this approach is that we only have to loop through the third molecule/atom and we only do this after we have determined that the first two are within the cut-off. Another (more intuitive) option would be to call Compute_ThreeBody_Energy within the move_* files (e.g., move_delete) after calling the other Compute_* subroutines (e.g., Compute_Molecular_Bond_Energy). From what I can tell, the downside to this approach is that we would be repeating the time consuming calculation of looping through N molecules, which is already being performed within Compute_Molecule_Nonbond_Inter_Energy.

In brief, do you think this approach would work? If you have concerns, what would be the simplest, fastest, or best way to implement a three-body potential in Cassandra?

Thank you

Eliseo,

Thank you for the response. Here are my thoughts regarding the two methods:

1) Yes, a neighbor list is definitely a good way to go. If you already have some code that would be useful for this, then I would certainly like to take a look at it. What would be the best way for you to share that with me? I guess I assumed that Cassandra already implemented a neighbor list, but is that not the case? From what I remember, neighbor lists are less common in Monte Carlo codes than Molecular Dynamics because of regrowth, insertion, and deletion moves. Is that why Cassandra does not use a neighbor list and instead stores the pair_nrg_array?

2) So essentially triad_nrg_array would be an NxNxN matrix? Do you have any concerns with regards to memory and speed if N >> 1000? As far as the implementation, are there any simplifications that we can make considering we only have one molecule type and that molecule consists of only a single atom (Helium)? Also, do we need to add two loops or just one? I thought we already have a double loop which goes over all pairs of interactions, so it seems like we would only need one additional loop. But perhaps I am misunderstanding your explanation. To address my confusion, can you briefly define what ispecies and nspecies represent? For example, a simulation with united-atom representations of 100 propane molecules and 50 butane molecules. Is nspecies 2, 100, 50, 150, (3*100), (4*50), or (3*100 + 4*50)? Or for the specific problem at hand, with a simulation of 1000 Helium atoms is nspecies 1 or 1000?

Thanks again

Rich